#ifndef MaTRIX_H
#define MaTRIX_H

#include <iostream>
#include "Vertex.h"

#define O_mult_matrix(res, arg1,arg2) typeof(arg1 * arg2) res = arg1 * arg2

namespace Y{

template <int M, int N>
class Matrix{

	public:

		float data[M][N];

	public:
		
		Matrix () { }
		Matrix ( float p[] ){
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < N; ++j){
					data[i][j] = p[j * M + i];
				}
			}
		}
		Matrix ( float p[M][N] ){
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < N; ++j){
					data[i][j] = p[i][j];
				}
			}
		}
		Matrix ( double p[M][N] ){
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < N; ++j){
					data[i][j] = p[i][j];
				}
			}
		}
		
		Matrix ( const Matrix<M,N>& t2 ){
			for ( int i = 0; i < M; ++i ){ for (int j = 0; j < N; ++j){ data[i][j] = t2(i,j); } }
		}
		Matrix<M,N>& operator= ( const Matrix<M,N>& t2 ){
			for ( int i = 0; i < M; ++i ){ for (int j = 0; j < N; ++j){ data[i][j] = t2(i,j); } }
			return *this;
		}

		template <int P>
		Matrix<M,P> operator* (const Matrix<N,P>& m){
	
			Matrix<M,P> nData;
			
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < P; ++j){
					nData(i,j) = 0.0;
					for (int k = 0; k < N; ++k){
						nData(i,j) += m(k,j) * this->get(i,k);
					}
				}
			}
			
			return nData;
		}
		
		
		template <int P>
		Matrix<M,P> operator* (float m[N][P]){
	
			Matrix<M,P> nData;
			
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < P; ++j){
					nData(i,j) = 0.0;
					for (int k = 0; k < N; ++k){
						nData(i,j) += m[k][j] * this->get(i,k);
					}
				}
			}
			
			return nData;
		}
				
		Matrix& operator*= (float f){ for ( int i = 0; i < M; ++i ){ for ( int j = 0; j < N; ++j ){ data[i][j] *= f; } } }
		
		inline float get (int i, int j){ return data[i][j]; }
		inline float& operator() (int i, int j){ return data[i][j]; }
		inline float operator() (int i, int j) const{ return data[i][j]; }
		
		void print (){

			std::cout << "\n";
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < N; ++j){
					std::cout << "| " << data[i][j] << " |";
				}
				std::cout << "\n";
			}

		}
		
		void makeIdentityMatrix() { for ( int i = 0; i < M; ++i ){ for ( int j = 0; j < N; ++j ){ if ( i == j ){ data[i][j] = 1.0; } else{ data[i][j] = 0.0; } } } }
		void loadIdentity();
		void translate(float dx, float dy, float dz);
		void translate(const Vertex3f&);
		void scale(float sx, float sy, float sz);
		void rotate_x(float deg);
		void rotate_y(float deg);
		void rotate_z(float deg);
		void negate_x();
		void negate_y();
		void negate_z();
};


		template <int M, int N, int P>
		static inline void mult (Matrix<M,P>& nData,   const Matrix<M,N>& m1 , const Matrix<N,P>& m){
			
			float buf[M][P] = {0.0};
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < P; ++j){
					for (int k = 0; k < N; ++k){
						buf[i][j] += m.data[k][j] * m1.data[i][k];
					}
				}
			}
			
			for ( int i = 0; i < M; ++i ){
				for ( int j = 0; j < P; ++j ){
					nData.data[i][j] = buf[i][j];
				}
			}

		}
		
		template <int M, int N, int P>
		static inline void mult (Matrix<M,P>& nData,   const Matrix<M,N>& m1 , float m[N][P]){
			
			float buf[M][P] = {0.0};
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < P; ++j){
					for (int k = 0; k < N; ++k){
						buf[i][j] += m[k][j] * m1.data[i][k];
					}
				}
			}
			
			for ( int i = 0; i < M; ++i ){
				for ( int j = 0; j < P; ++j ){
					nData.data[i][j] = buf[i][j];
				}
			}

		}
		
		template <int M, int N, int P>
		static inline void mult (Matrix<M,P>& nData,   float m1[M][N], const Matrix<N,P>& m){
			
			float buf[M][P] = {0.0};
			for (int i = 0; i < M; ++i){
				for (int j = 0; j < P; ++j){
					for (int k = 0; k < N; ++k){
						buf[i][j] += m.data[k][j] * m1[i][k];
					}
				}
			}
			
			for ( int i = 0; i < M; ++i ){
				for ( int j = 0; j < P; ++j ){
					nData.data[i][j] = buf[i][j];
				}
			}

		}
		
		#include "MatrixTransformations.h"
}
#endif